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SECTION  I 
INTRODUCTION 


An  example  of  an  important  application  of  three-dimensional  Euler  equation 
solutions  on  dynamic  grids  is  that  of  transonic  flow  about  non-stationary 
stores.  Unsteady  aerodynamic  information  is  extremely  difficult  and  costly  to 
obtain  in  wind  tunnel  experiments,  particularly  in  the  transonic  regime. 
Consequently,  it  is  important  to  know  when  steady-state,  quasi-steady,  or  true 
unsteady  aerodynamic  information  is  required.  Computational  fluid  dynamics 
(CFO)  has  matured  to  the  point  where  it  can  be  used  to  assist  in  deciding  the 
conditions  for  which  steady  or  unsteady  information  is  needed,  as  well  as  being 
itself  a  source  for  providing  such  information. 

To  compute  the  flow  of  a  store  released  from  an  aircraft,  it  is  desirable 
to  solve  the  unsteady  Euler  equations  on  a  grid  that  moves  with  the  store  along 
its  trajectory.  The  purpose  of  this  paper  is  to  solve  the  three-dimensional 
unsteady  Euler  equations  on  a  time-dependent  grid.  The  computations  presented 
here  are  for  bodies  whose  motion  is  prescribed.  As  the  solution  advances  in 
time,  body  motion  could  be  determined  from  the  Euler  equations  by  using  force 
and  moment  coefficients  obtained  from  the  Euler  solution  in  the  dynamic  equa¬ 
tions  of  motion  for  the  body  to  determine  the  trajectory  of  the  body.  The 
objective  of  this  report  is  to  present  and  verify  dynamic-grid  Euler  equations 
computations . 

The  Euler  equations  solution  method  used  is  similar  to  that  presented  by 
Janus,  Reference  1.  It  is  a  finite  volume,  flux-vector  split,  second-order 
explicit  scheme,  that  can  use  local  time  stepping  with  CFL  <2  for  steady-state 
computations,  Reference  2.  The  flux  vectors  at  cell  faces  are  determined 
differently  in  the  present  method.  In  References  1  and  2  the  flux  vectors  at 
cell  faces  depend  on  a  cell  face  eigenvalue  that  was  determined  from  informa¬ 
tion  resulting  from  averaging  dependent  variables  on  either  side  of  the  cell 
face.  Due  to  the  splitting,  each  flux  vector  is  the  sum  of  three  subvectors, 
each  of  which  has  an  eigenvalue  as  a  coefficient.  The  cell  face  eigenvalue, 
computed  as  described  above,  was  used  to  compute  a  flux  vector;  however,  the 
remaining  portion  of  the  flux  vector  was  computed  with  information  from  the 


appropriate  direction  as  dictated  by  the  sign  of  the  eigenvalue.  In  the 
present  method,  each  subvector  of  a  flux  vector  is  computed  on  either  side  of  a 
cell  face  by  extrapolating  the  dependent  variables  from  either  side  of  a  cell 
face.  An  eigenvalue  on  either  side  of  a  cell  face  is  computed  in  like  manner. 
The  flux  vector  at  a  cell  face  is  then  formed  by  summing  the  subvectors  on 
either  side  of  a  cell  face  by  appropiately  taking  into  account  the  sign  of  the 
eigenvalues.  This  method  requires  less  storage  and  executes  more  quickly  than 
the  previous.  It  also  totally  eliminates  upstream  propagation  of  information 
in  supersonic  regions.  Experience  with  this  technique  has  shown  that  the 
ringing  downstream  of  shocks  characteristic  of  the  previous  method  is  greatly 
reduced.  '  disadvantage  of  the  method  is  that  in  practice  it  frequently 
requires  a  slightly  lower  CFL  number  for  stability  than  the  previous  method. 

Characteristic  variable  boundary  conditions  (Reference  1)  are  used  on  all 
boundaries  with  the  exception  of  the  reflection  plate  used  in  the  moving  store 
calculations.  On  the  reflection  plate,  reflection  plane  boundary  conditions 
are  used. 


SECTION  II 
ALGORITHM 

1.  BASIC  EQUATIONS 

The  conservation  law  form  of  the  Euler  equations  in  Cartesian  coordinates 
x,  y,  and  z  is 

(.) 

where 

q  •  [p,  pu,  pv,  pw,  e]T 

f  ■  [pu,  pu2  ♦  p,  puv,  puw,  u(e+p)]^ 

g  -  [pv,  puv,  pv2  ♦  p,  pvw ,  v(e+p)]T 

h  -  [pw,  puw,  pvw,  pw2  ♦  p,  w(e+p ) ]T 

P  "  (Y  ~  1)[e  -  p(u^  +  v2  ♦  w~)] 

Using  time-dependent  curvilinear  coordinates  defined  by 

£  "  £(x,  y.  z,  t) 

n  *  n(x,  y,  2,  t ) 


C  -  c(x,  y,  2,  t) 


Equation  (1)  may  be  transformed  to  (Reference  3) 


IQ  ♦  K  ♦  i£  ♦  ii!  .  o 

3t  3£  3n  3; 


where 


Q  -  J[p.  pu,  ov,  pw,  e] 


F  -  J[pU,  pull  ♦  £xp,  pvU  *  £yp,  pwl)  ♦  £„p,  U(e*p)  -  t.p]^ 


0  -  JCpV.  puV  ♦  nxp,  pvV  +  n  P.  pwV  ♦  nzp,  V(e+p)  -  ntp]1 


H  -  J[pW,  puW  ♦  cxp,  pvW  ♦  cyp,  pwW  ♦  C2P.  W(e+p)  -  c^p] 


J  -  x^y^z  ~  Zr)yz)  '  y£(xnzc  "  2nxc)  *  Zt(Xny;  "  yr<-' 


"  J_1{ynrc'*nyC) 

nx  ■  J'1(zzyc~ys2^ 


*  J"  (y«  Wn* 


“  J"  (xnyc'ynxC) 


\  *  J"  (Vc‘yeV 


cz  "  J"  (xCyn"yCxn) 


Ky  m  J’I(Vc'Vc) 


by  "  3  ^XtZC_Z5XC^ 


cy  -  j'  (xnzrznxe) 


“  _XT^x”yT^y  ZT^Z 

\  ’  -XTnx-yTny'*Tnz 

*t  "  ‘XTVyTVZTl* 


n  -  cxu  +  ?yv  +  $zw  +  £t 
v=nu+nv+nw+n 

x  y  z  t 

W-CU  +  CV  +  C.w  +  C,. 

x  y  z  c 


A  finite  volume  discretization  of  Equation  (2)  balances  the  increase  of 
the  conserved  quantity  in  a  computational  cell,  or  volume,  with  the  flux  of  the 
quantity  through  the  surface  of  the  cell.  Figure  1  depicts  a  portion  of  the 
computational  domain  with  a  typical  cell  labeled.  Assuming  the  dependent 
variables  are  constant  in  the  interior  of  cell  i,  j,  k,  and  that  the  flux 
vectors  F,  G,  and  H  are  constant  over  the  constant  £  ,  constant  n  and  constant  K 
surfaces  of  the  cell,  respectively,  an  explicit  discretization  of  Equation  (2) 
is 
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This  can  be  rewritten  as 


AQ.  6F  JO.  AH, 

.  _1  .  _i  .  _Ji  .  o 

At  AC  An  A( 


(«) 


Note  that  any  algorithm  written  in  this  form  is  conservative  as  long  as  the 
fluxes  are  sing  1  e-  val  ued  functions  of  the  cell  faces. 


v 


f  =  Constant  Plane 

Figure  1.  Computational  Grid 


2.  FLUX  SPLITTING 


Hyperbolic  partial  differential  equations,  such  as  the  Euler  equations,  are 
characterized  by  the  existence  of  a  limited  domain  of  dependence.  The  solution 
at  a  point  does  not  depend  on  every  other  point  in  the  field;  this  means  that 
information  travels  only  in  certain  characteristic  directions.  Numerical 
schemes  intended  to  solve  hyperbolic  equations  are  usually  enhanced  by  insuring 
that  the  numerical  method  propagates  information  in  the  direction  specified  by 
the  partial  differential  equation.  This  can  be  done  by  using  an  upwind  method, 
or  one  in  which  the  difference  operator  is  taken  in  the  direction  from  which 
the  information  should  come.  Stability  properties  are  often  improved  by 
upwinding,  and  it  is  usually  unnecessary  to  add  smoothing  terms  or  artificial 
viscosity  to  an  upwind  method. 


The  three-dimensional  Euler  equations.  Equation  (2),  are  a  hyperbolic 
system  of  five  equations  and  hence  have  five  characteristic  velocities  in  each 
of  the  three  spatial  directions.  These  character istic  velocities  are  deter¬ 
mined  from  the  quasi  linear  form  of  Equation  (2), 


3Q 

3t 


♦ 


(5) 


where  the  matrices  A,  B,  and  C  are  given  by 


A  - 


3F 

55* 


B 


3G 
"  3Q' 


C 


3H 

3Q 


The  eigenvalues  of  A  are  the  characteristic  velocities  in  the  £  direction;  the 
eigenvalues  of  B  are  the  character istic  velocities  in  the  v  direction;  and 
similarly  the  eigenvalues  of  C  are  the  characteristic  velocities  in  the  T 
direction . 


Since  F,  G,  and  H  are  identical  except  that  where  £  appears  in  F,  n  will 
appear  in  G  and  C  will  appear  in  H,  extra  writing  can  be  avoided  by  letting  K 
represent  either  F,  G,  or  H,  and  k  represent  either  £,  >j ,  or  f  . 

Then  define 


which  corresponds  to  A,  B,  or  C  depending  on  the  meaning  of  K.  The  eigenvalues 
of  the  matrix  K  are 


Xk  “  Xk  *  Xk  “  kxu  +  V  +  k*W  +  kt  "  bk 

Xk  “  6k  +  c|Vk|  "  6k  +  C(kx  +  yl  +  (?) 

Xk  "  6k  “  C,7k| 


where  c  is  the  speed  of  sound. 


It  is  possible  to  split  the  flux  vector  K  into  three  parts,  one  corre¬ 
sponding  to  each  of  the  distinct  eigenvalues  of  K  given  above.  (For  the 
details  of  this  splitting,  see  Reference  1).  The  flux  vector  K  is  then  written 
as 


K 

- 

U  +  XkK5 

(8) 

where 

’  ' 
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P 
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pu+pci<x 

pv 

V  J 

*4  "  2Y 

pv+pcky 
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f(u2*v2*v2) 

e+p+pce 

*5  “  2Y  I  pV_pckv 


pw-pck 


p»p-pc0^. 


,.2.2.  2.2 
(k  +k  +k  ) 
x  y  z 


ku+kv+kw. 
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The  sign  of  A<in  Equation  (8)  determines  from  which  direction  information 
should  be  used  to  determine  the  corresponding  portion  of  the  flux,  K^. 

3.  NUMERICAL  IMPLEMENTATION 


For  the  finite  volume  discretization,  Equation  (4),  the  dependent  variables 
are  stored  at  cell  centers,  but  the  fluxes  are  required  at  cell  faces.  There¬ 
fore,  some  type  of  interpolation  or  extrapolation  must  be  used  to  determine  the 
fluxes.  The  signs  of  the  eigenvalues  at  the  faces  should  determine  the  direc¬ 
tion  of  the  extrapolation  to  be  used,  but  note  that  there  is  still  some 
ambiguity,  since  the  eigenvalues  themselves  may  be  calculated  in  several 
different  ways  at  cell  faces.  Whitfield  and  Janus  (Reference  2)  computed  the 
eigenvalues  at  a  face  from  the  average  of  dependent  variables  in  the  two  cells 
sharing  the  face.  The  present  scheme  uses  a  different  approach.  The  scheme 
reported  here  calculates  a  set  of  left  eigenval ues,  a* (Q*-) ,  and  left  split 


fluxes,  K*{Ql)  from  dependent  variables  extrapolated  from  cell  centers  left  of 
the  face,  QL,  and  it  also  calculates  a  set  of  right  eigenval ues, (Qr),  and 
right  split  fluxes,  Kj,(Qr),  from  dependent  variables  extrapolated  from  cell 
centers  right  of  the  face,  QR.  The  flux  at  the  face  is  then  set  to 

K(QL.QR>  -  I  (QL)  ♦  |aJ(QL)|)K  (QL) 

fl.4.5  1 


♦  j(A*(QR)  -  |aJ(QR)|)Kji(QR)] 


Note  that  if  a  left  and  right  eigenvalue  have  different  signs,  then  the  corre¬ 
sponding  split  flux  may  be  either  summed  from  both  sides,  if  the  sign  changes 
from  positive  to  negative,  or  from  neither  side,  if  the  sign  changes  from 
negative  to  positive. 


Similar  extrapolations  are  used  in  the  ij  and  f  directions.  Taking  A(sAnsAf 
the  predictor  is 

5i,j,k  "  Qi,J,k  "  ATl,j,kt4*i  *  a6j  (13) 


where 


**t  “  F<«L  ,  •  QR  ,  )  -  r( ql  .  q"  .  ) 


and  6Gf.  and  6H  are  evaluated  similarly. 
J  K 


The  corrector  is  written  as 


Qn+1  ■  Qn  -  — LbljL^kiF  *40  ♦  4H  ] 
l.J.k  yl,J,k  2  Lorl  0V*J  onkJ 


where 


1»^C4?1  ♦  40 j  ♦  4^] 


*  *■  *1 
SF  -  F(QL  q"  )  -  r(QL 

i*j.J.k,  i*j,J,k  1 


L  *QR  ) 

l"2»J»k.  i“j>»  J  *k 


and  6Gj  and  6H^  similarly  are  flux  differences  using  two-point  dependent 
variable  extrapolations.  6^.,  68.,  and  represent  differences  of  fluxes 
evaluated  using  a  one-point  extrapolation  of  the  predicted  variables. 

The  time  step  Ar-jj^  is  given  by 


lijik  1  »J  ,k  1  tj  »k  i  fj  »k  i,J,k  i»J»k 


for  k  =  t,i)  ,  and  f .  When  time  accuracy  is  not  required  (steady  state)  the 
maximum  time  step  for  each  individual  cell  is  used  where  CFL  <  2.  This  is  done 
to  accelerate  convergence  by  enabling  every  cell  to  advance  in  time  as  fast  as 
possible.  The  edge  cells  (those  in  the  domain  adjacent  to  a  boundary)  are  only 
first-order  accurate  and  stable  for  a  CFL  is  1,  consequently  the  time  steps  for 
these  cells  are  reduced  accordingly.  In  situations  where  a  time  accurate 
solution  is  desired,  such  as  in  some  dynamic  grid  applications,  the  minimum  of 
all  the  local  time  steps  is  selected  for  use  at  all  cells.  This  enables  the 
solution  to  advance  at  the  same  time  in  every  cell,  although  restricting  the 
time  step  to  that  of  the  slowest  cell.  This  method  of  selecting  the  time  step 
is  referred  to  as  minimum  time  stepping. 

4.  BOUNDARY  CONDITIONS 

Consistent  with  the  impetus  leading  to  the  flux  split  form  of  the  Euler 
equations,  i.e.,  honoring  propagation  directions,  a  set  of  boundary  conditions 
has  been  developed,  referred  to  as  the  characteristic  variable  boundary  condi¬ 
tions  (CVBC) ,  Reference  2.  The  CVBC  presented  in  Reference  2  are  formulated 
for  stationary  grids.  Some  minor  adjustments  to  these  equations  are  required 
to  establish  a  valid  set  of  CVBC  for  the  dynamic  grid  case. 

The  characteristic  variables  are  constructed.  Reference  1,  using  the  matrix 
of  eigenvectors  of  the  flux  Jacobian  and  the  non-conservative  dependent  vari¬ 
able  vector.  Since  in  the  dynamic  grid  derivation  neither  the  eigenvectors  or 
dependent  variables  are  altered.  Reference  1,  the  characteristic  variables 
remain  unchanged.  Consequently  the  equations  corresponding  to  the  cases 
described  in  Reference  2  are  all  valid  except  for  the  case  of  the  impermeable 
surface.  This  case  will  be  reformulated  for  the  dynamic  grid  application. 

The  impermeable  surface  case  is  characterized  by  the  first  three  eigen¬ 
values  equal  to  zero,  the  fourth  positive,  and  the  fifth  negative.  For  this 


case  no  flow  through  the  surface  must  be  specified,  replacing  one  of  the 
characteristic  variable  relation  equations  yielding 
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where  the  r  subscripts  refer  to  a  reference  value,  which  is  selected  as  the 
center  of  the  first  cell  from  the  boundary,  and  the  b  subscripts  refer  to 
boundary  values.  Solving  Equations  17d,e  for  Ap  yields 


4p  -  pb  -  pr  .  ;  p0c0(kxur  ♦  kyvr  ♦  k2wp  ♦  ktfb>  (18a) 

The  other  relations  can  be  obtained  by  solving  Equations  (17a,b,c,d)  simulta¬ 
neously  yielding 
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The  minus  and  plus  signs  in  Equations  (18a)  correspond  to  the  location  of  the 
point  r.  ff  the  point  r  is  in  the  positive  k  direction  from  the  boundary  then 
the  minus  sign  is  used,  and  if  it  is  in  the  negative  direction  then  the  plus 


SECTION  III 
RESULTS 

Dynamic-grid  computations  were  performed  for  the  following  configurations: 
(a)  plunging  and  oscillating  airfoils,  (b)  plunging  wings,  (c)  horizontally 
launched  missile  (or  store),  and  (d)  vertically  launched  store  from  a  reflec¬ 
tion  plate.  The  plunging  airfoil  and  wing  cases  were  run  for  the  special  case 
of  constant  velocity  plunge  in  order  to  check  the  dynamic-grid  code.  A  con¬ 
stant  velocity  plunge  corresponds  to  a  stationary  solution  at  some  particular 
angle  of  attack.  The  constant  velocity  plunge  solutions  were  compared  with 
steady-grid  angle  of  attack  solutions  to  insure  that  both  solutions  were  iden¬ 
tical,  which  they  were.  The  same  was  done  for  the  horizontally  launched 
missile.  In  one  case  the  missile  was  launched  into  still  air  and  then  the 
solution  was  compared  with  a  solution  for  a  stationary  missile  with  the  air 
flowing  past.  These  results  were  reported  in  Reference  1.  The  same  numerical 
experiment  was  also  performed  for  a  store  in  constant  velocity  vertical  plunge 
and  a  stationary  store  at  angle  of  attack.  Again  the  solutions  were  identical 

The  cases  selected  are  the  solution  for  an  oscillating  airfoil  and  the 
solution  for  a  vertically  launched  store  from  a  reflection  plate.  The  oscil¬ 
lating  airfoil  was  investigated  because  of  the  need  to  determine  whether  the 
phase  was  being  computed  correctly.  Experimental  data  are  available  for  an 
NACA  0012  airfoil.  Reference  5,  consequently  this  airfoil  was  selected.  The 
experiment  consisted  of  the  harmonic  pitch  oscillation  of  a  NACA  0012  airfoil 
about  the  quarter  chord  point.  Figure  1  is  the  result  obtained  from  the  har¬ 
monic  oscillation  using  a  reduced  frequency  of  0.2156  and  oscillating  *  2.51° 
about  a  mean  alpha  of  0.016°  in  a  uniform  Mach  0.755  flow.  In  an  attempt  to 
isolate  the  phase  shift  of  the  lift  curve,  the  normalized  lift  curves  of  the 
numerical  solution  and  experimental  data.  Reference  3,  are  shown  in  Figure  2. 

It  is  evident  from  this  figure  that  the  correct  phase  is  calculated.  The  lift 
curves  could  not  be  directly  compared  because  the  numerical  solution  was 
obtained  using  a  rather  coarse  grid,  consequently  the  resolution  of  the  values 
of  lift  and  drag  is  not  adequate. 
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Figure  2.  Phase  Shift  (computed  and  experimental) 


Finer  grid  steady-state  solutions  were  obtained  and  compared  to  previous  lift 
and  drag  Euler  solutions  to  make  sure  the  present  code  produces  acceptable 
results  for  steady  state.  For  example,  for  Mot.  =  0.8  and  a  =  1.25°  the  lift  and 
drag  coefficients  are  C|_  =  0.3796  and  Cq  =  0.0211,  whereas  Anderson,  Thomas, 
and  Van  Leer,  Reference  6,  give  CL  =  0.3740  and  Cq  =  0.0232. 

The  second  case  presented  here  is  that  of  a  store  launched  vertically  from 
a  reflection  plate.  There  appear  to  be  no  experimental  data  available  for 
stores  moving  in  the  vicinity  of  wings  or  other  stores.  There  are  steady-state 
data  available  for  a  store  located  at  various  positions  relative  to  a  reflec¬ 
tion  plate.  Reference  7.  A  grid  with  60  points  in  the  streamwise  direction,  10 
points  in  the  radial  direction,  and  30  points  in  the  circumferential  direction 
was  used  to  obtain  comparisons  with  the  steady  experimental  data  for  M«  =1.41 
and  a  zero  store  angle  of  attack  with  respect  to  the  reflection  plate.  The 
model  geometry  and  grid  are  illustrated  in  Figure  3.  A  plane  of  symmetry 


perpendicular  to  the  reflection  plate  and  bisecting  the  store  was  taken  advan¬ 
tage  of  to  half  the  size  of  the  problem.  A  comparison  is  presented  in  Figure  4 
of  computed  and  measured  store  surface  pressure  distributions  for  a  store  nose 
to  reflection  plate  separation  distance  of  10.07  percent  of  the  store  length. 
The  dip  in  the  Cp  curve  near  the  nose  is  due  to  the  reflection  of  the  shock 
emanating  from  the  store  nose.  Although  viscous  effects  are  not  accounted  for 
on  either  the  store  or  reflection  plate  in  the  present  computations,  good 
qualitative  agreement  is  obtained  between  the  computed  and  measured  steady- 
state  data. 


Figure  3  (a).  Perspective  View  of  Test  Setup 
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Figure  4.  Computed  and  Experimental  Store  Surface  Pressures 


To  help  assess  the  importance  of  unsteady  aerodynamics  during  store  separa¬ 
tion,  dynamic-grid  Euler  solutions  were  obtained  for  the  store  discussed  above 
when  launched  from  the  10.07  percent  location  in  a  direction  perpendicular  to 
and  away  from  the  plate.  The  procedure  used  for  this  calculation  was  to  first 
calculate  the  initial  steady-state  flow  field  with  the  store  at  the  10.07 
percent  location.  The  plunge  velocity  of  the  store  away  from  the  wall  was  then 
suddenly  set  to  a  predetermined  value  and  held  constant  thereafter.  Present- 
day  store  separation  analyses  typically  use  steady  aerodynamics  with  a  correc¬ 
tion  to  the  store  angle  of  attack  to  account  for  the  unsteady  store  motion. 

The  geometric  angle  of  attack  of  the  store,  i.e.  the  pitch  angle  e  in  Figure 
3(a),  is  incremented  by  the  effective  angle  of  attack  due  to  the  instantaneous 


store  plunge  velocity.  In  other  words,  if  the  free  stream  velocity  is  tL  end 
if  the  store  has  a  plunge  velocity  of  v,  then  the  aerodynamic  forces  for  this 
case  are  usually  obtained  from  steady  aerodynamics  with  c  incremented  by  an 
amount  given  by 

au  -  tan-1  (19) 

To  evaluate  the  validity  of  this  approximation,  a  set  of  steady  computations 
were  performed  for  a  store  with  e  =3  degrees  and  compared  with  two  unsteady 
dynamic-grid  computations.  For  these  calculations,  a  60  by  10  by  15  grid  was 
used.  The  first  dynamic-grid  computation  had e  =0  degrees  and  a  downward  veloc¬ 
ity  such  that  the  effective  angle  of  attack  due  to  the  velocity,  uu  ,  was  3 
degrees.  The  second  dynamic-grid  computation  had  e  =2  degrees  and  an  effective 
angle  of  attack  increment  due  to  velocity  of  one  degree.  As  explained  earlier, 
the  aerodynamics  for  both  of  these  unsteady  cases  would  normally  be  modelled  in 
a  store  separation  analysis  by  a  steady-state  condition  with  e  =3  degrees.  A 
comparison  of  the  store  pressure  distribution  along  the  0=0  plane  (see  Figure 
3)  for  the  steady  and  two  unsteady  cases  is  presented  in  Figure  5.  The  reflec¬ 
tion  plate  pressure  distribution  for  the  same  conditions  is  given  in  Figure  6. 
From  the  figures  it  is  apparent  that  a  noticeable  difference  does  exist  between 
the  steady  and  unsteady  cases  particularly  during  the  initial  moments  of  the 
store  motion.  Note  that  the  pressures  on  the  store  ahead  of  the  reflected 
shock  for  the  au  =  i  degree  case  quickly  become  indistinguishable  for  the  steady 
results,  but  the  au  =3  case  shows  a  larger  difference.  The  greatest  difference 
in  results  is  due  to  the  position  of  the  reflected  shock.  A  comparison  of 
normal  force  and  pitching  moment  coefficients  is  given  in  Figure  7.  Figure 
7(a)  shows  a  spike  in  the  normal  force  coefficient  at  the  first  time  step  which 
is  due  to  apparent  mass  effects  from  the  impulsive  change  in  velocity.  This  is 
a  transient  effect  which  decays  quickly.  After  the  decay  of  the  large  tran¬ 
sient  due  to  store  acceleration,  a  difference  between  the  steady  and  unsteady 
cases  still  exists.  This  is  due  at  least  in  part  to  the  lag  in  the  reflected 
shock's  position  for  the  unsteady  case. 
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Figure  5.  Comparison  of  Store  Surface  Pressure  Along  6=0  for  Unsteady 

Store  Motion  and  For  Steady  Store  at  Equivalent  Angle  of  Attack 
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Figure  7(a).  Store  Normal  Force  Coefficient 
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Figure  7  (b) .  Coefficient  of  Pitching  Moment  About  59.8  Percent 


These  results  illustrate  the  difference  in  aerodynamic  forces  for  a  pres 
cribed  store  trajectory.  It  is  possible  to  determine  the  store  trajectory 
using  forces  and  moments  obtained  from  the  Euler  equations,  as  Mastin  did. 
(Reference  8)  To  determine  the  significance  of  true  unsteady  aerodynamics 
to  store  separation,  the  rigid  body  equations  of  motion  should  be  solved 
using  both  the  unsteady  and  the  steady  aerodynamics  so  that  the  difference 
in  the  store  trajectory  can  be  seen. 


SECTION  IV 
CONCLUDING  REMARKS 


A  method  was  presented  for  solving  the  three-dimensional  unsteady  Euler 
equations  on  dynamic  grids  based  on  flux  vector  splitting.  The  equations  were 
cast  in  curvilinear  coordinates  and  a  finite  volume  discretization  was  used  for 
handling  arbitrary  geometries.  The  discretized  equations  were  solved  using  an 
explicit  upwind  second-order  predictor-corrector  scheme  that  is  stable  for  a 
CFL  of  2.  Characteristic  variable  boundary  conditions  were  developed  and  used 
for  unsteady  impermeable  surfaces  and  for  the  far-field  boundary.  Dynamic-grid 
results  were  presented  for  an  oscillating  airfoil  and  for  a  store  separating 
from  a  reflection  plate. 

It  was  determined  that  for  the  cases  considered  of  stores  separating  from  a 
reflection  plate  that  the  unsteady  aerodynamic  forces  on  the  store  are  signifi¬ 
cantly  different  from  forces  obtained  by  steady-state  aerodynamics  with  the 
body  inclination  angle  changed  to  account  for  plunge  velocity. 
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